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Abstract 

In  this  thesis  we  consider  problems  for  which  the  boundary  is  not  known  before 
the  problem  is  solved  and  must  be  determined  as  part  of  the  solution.  We  consider 
a  time  dependent  problem  which  results  in  a  moving  boundary.  We  look  at  the 
heat  conduction/diffusion  equation  in  one  and  two  spatial  dimensions.  We  introduce 
the  fundamental  solution  or  Green’s  function  and  use  Green’s  Theorem  to  yield  a 
Volterra.  boundary  integral  equation  which  involves  an  unknown  function  on  the  mov¬ 
ing  boundary.  We  take  the  limit  of  our  integral  expression  to  the  moving  boundary 
to  obtain  a  nonlinear  system  of  integral  equations  for  the  location  of  the  boundary 
and  the  unknown  function.  We  use  the  boundary  element  method  to  obtain  a  so¬ 
lution  to  this  system  of  integral  equations.  This  solution  is  then  substituted  back 
into  the  original  Volterra  equation  to  obtain  the  solution  of  our  original  problem. 
Graphical  results  for  the  two  dimensional  problem  are  presented. 
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THE  BOUNDARY  ELEMENT  METHOD  APPLIED 
TO  THE  TWO  DIMENSIONAL  STEFAN 
MOVING  BOUNDARY 
PROBLEM 


I.  Introduction 


1.1  Back  round 

Typically,  the  goal  in  solving  differential  equations  is  to  find  a  solution  of  the 
equation  on  some  domain  0  given  the  value  of  the  solution  or  the  value  of  the  normal 
derivative  of  the  solution  on  the  boundary  T  of  the  domain.  But,  there  is  a  class 
of  boundary  value  problems  for  which  the  boundary  is  not  known  in  advance  and 
must  be  determined  as  part  of  the  solution.  These  problems  are  known  as  unknown 
boundary  problems.  If  the  differential  equation  is  steady  state  (no  time  derivatives) 
then  the  boundary  will  be  unknown  but  will  remain  fixed.  These  are  known  as 
free  boundary  problems.  If  the  equation  is  evolutionary  then  the  boundary  will 
move  as  a  function  of  time  and  the  domain  on  which  we  are  solving  the  differential 
equation  will  be  time-dependent.  These  are  known  as  moving  boundary  problems. 
This  thesis  will  deal  with  a  two-dimensional  moving  boundary  problem.  In  order 
to  solve  a  problem  with  an  unknown  boundary,  we  must  impose  extra,  conditions 
on  the  solution's  behavior  at  the  boundary  in  addition  to  the  typical  boundary  and 
initial  conditions.  These  extra  conditions  are  often  drawn  from  physical  concerns 
such  as  the  conservation  of  energy  or  the  physical  spatial  constraints  imposed  on  a 
real-world  problem. 

Unknown  boundary  problems  arise  frequently  in  the  physical  world.  They  can 
be  found  in  the  study  of  the  solidification  of  metals,  carbon  diffusion  in  steelmaking. 
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ablation  of  materials  using  a  laser  beam,  welding  of  two  metals,  corrosion  and  oxi¬ 
dation  of  metals,  oxygen  diffusion  in  biological  tissues,  the  melting  and  freezing  of 
ice,  thermal  switching  of  glasses,  stellar  evolution,  diffusive  chemical  reactions,  elec¬ 
trochemical  machining,  water  seepage  through  darns,  impulse  and  optimal  control, 
and  crystal  growth.  This  thesis  will  look  at  an  example  involving  the  melting  and 
freezing  of  ice. 

When  considering  the  melting  of  ice  it  is  possible  to  have  a  single  phase  problem 
or  a  two  phase  problem.  In  the  single  phase,  the  ice  is  initially  at  it’s  melting 
temperature  V  =  0°O  and  only  the  heat  of  fusion  need  be  transferred  in  order  to 
melt  it.  In  the  two  phase  problem,  the  ice  is  initially  at  some  temperature  below  it's 
melting  point  and  we  must  transfer  heat  energy  in  order  lo  bring  it  to  the  melting 
point  and  then  supply  additional  heat  energy  in  order  to  satisfy  the  heat  of  fusion. 
VVe  will  consider  the  single  phase  problem  in  this  thesis. 

1.2  Coals 

The  goals  of  this  thesis  are: 

1 .  To  model  phase  change  problems  in  one  and  two  dimensions. 

2.  To  apply  the  boundary  element  technique  to  a  two  dimensional  melting  ice 
problem. 

3.  To  use  computer  graphics  as  a  tool  to  display  the  numerical  results  in  a.  visual 
format. 

1.3  Scope 

We  have  limited  ourselves  to  two  dimensions  and  to  a  circular  initial  domain 
in  order  to  have  an  analytic  solution  available  for  comparison  with  the  numerical 
calculation  results.  Isotropic  materials  are  considered  to  simplify  the  heat  energy 
exchange.  Relatively  short  time  intervals  are  considered  to  avoid  excessive  calcula¬ 
tion  time.  Grid  refinement  is  kept  as  minimal  as  possible,  also  to  avoid  excessive 
computation  time.  We  neglect  changes  in  the  density  under  phase  transitions  in 
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order  to  avoid  dealing  with  convection  currents.  Convection  currents  will  occur  if 
the  density  is  allowed  to  change  even  if  the  liquid  phase  is  incompressible  and  its 
thermal  expansion  can  be  neglected  (3:227).  It  should  be  noted  that  our  model  can 
be  directly  applied  to  diffusion  problems. 

1.4  Support 

We  required  the  use  of  a  number  of  computer  systems  and  software  packages 
in  order  to  complete  this  thesis.  The  main  thrust  of  the  computation  was  borne  by 
the  A  FIT  ELXSl  computer.  The  large  number  of  test  runs  and  the  length  of  the 
runs  made  the  use  of  this  computer  a  necessity.  Shorter  length  computations  were 
performed  on  a  VaxStatiou  II/GPX  belonging  to  WL/ AARM.  Text  processing  was 
borne  by  AFIT's  Scientific  Support  Computer  VAX  11/785  and  also  by  a  WL/AARM 
VAX  11/785.  The  text  processing  was  accomplished  with  the  LNTj^X  Document 
Preparation  System. 

The  programming  code  was  written  in  the  Fortran  programming  language. 
This  language  ,vas  chosen  because  it  is  a  standard  programming  language  for  which 
optimized  compilers  are  rea  ;ily  available.  All  of  the  code  is  compliant  with  the  ANSI 
X-3.9  1978  (Fortran-77)  standard.  The  support  command  language  routines  were 
written  in  the  Csh  shell  programming  language.  Various  Unix  utilities,  such  as  awk 
and  sort,  were  used  in  the  shell  scripts. 

The  graphics  output  was  generated  by  subroutine  calls  to  the  metalib  graph¬ 
ics  package  on  the  Scientific  Support  Computer  and  by  the  CA-DISSPLA  graphics 
package  on  the  WL  VAXen.  The  metalib  software  is  a  local  adaptation  of  a  package 
that  was  developed  at  the  Air  Force  Weapons  Laboratory  at  Kirtland  AFB.  Local 
enhancements  were  made  by  Lt  Col  James  Lupo  of  AFIT/ENP.  The  A  FIT  Sun  work¬ 
station  network  was  utilized  in  running  the  Mathematica  software  for  calculations 
involving  Bessel  functions. 
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1.5  Mathematical  History  of  the  Problem 

Moving  boundary  problems  are  often  referred  to  as  Stefan  problems  due  to  the 
work  done  by  J.  Stefan,  in  1889,  on  the  freezing  of  the  ground  (21:173-484)  and  the 
melting  of  ice  layers  (22:965-983). 

Actually,  the  first  work  in  the  moving  boundary  area  was  performed  by  G. 
Lame  and  B.P.  Clapeyron  in  1831.  They  tried  to  calculate  the  thickness  of  a  solid 
crust  on  liquid  cooling  in  the  half  space  ,r  >  0  with  a  constant  temperature  on  the 
plane  x  =  0.  They  discovered  that  the  thickness  is  proportional  to  the  square  root 
of  time  but  did  not  arrive  at  the  constant  of  proportionality  (13:25 — 256). 

There  are  few  analytic  solutions  to  moving  boundary  problems  in  closed  form. 
Generally,  they  are  for  a  one-dimensional  geometry  on  an  infinite  or  semi-infinite 
domain  with  simple  boundary  and  initial  conditions  and  having  constant  thermal 
properties.  The  solutions  are  usually  functions  of  the  single  variable  .r/ \ft  and  are 
called  similarity  solutions.  Both  of  the  problems  that  Stefan  worked  on  possess 
similarity  solutions  of  the  form  u  =  2 <x\/i.  The  value  of  o  can  be  determined  from 
transcendental  equations  (17:2).  The  motive  behind  similarity  solutions  is  to  reduce 
the  number  of  independent  variables  by  taking  an  algebraic  combination  of  them 
(1:63-75). 

M.  Brillouin  reduced  the  solution  of  the  Stefan  problem  to  a  system  of  nonlinear 
integro-differential  equations  in  1929  (2:285-308).  However,  he  did  not  try  to  solve 
the  system  as  he  thought  it  would  be  very  difficult  to  do  so. 

In  1931,  L.S.  Leibenzon  derived  an  approximation  for  the  solution  of  the  Stefan 
problem  for  many  cases  (14:435-439).  lie  replaced  the  true  temperature  distribution 
by  a  quasi-stationary  solution.  This  solution  obeyed  the  Laplace  equation  in  space 
in  a  domain  that  had  a  moving  boundary  that  corresponded  with  the  solution  of 
Stefan. 
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Another  analytic  solution  technique  used  in  moving  boundary  problems  is 
transforming  the  coordinates  so  the  moving  boundary  becomes  fixed  in  the  new 
coordinate  system.  The  concept  of  conformal  transformations  is  used  in  the  hodo- 
graph  method  which  finds  frequent  application  in  the  field  of  fluid  mechanics.  The 
hodograph  method  is  presented  in  (4:288-293)  and  (15). 

Numerical  methods  applied  to  unknown  boundary  problems  have  an  extensive 
repertoire.  Since  the  solution  of  the  unknown  boundary  problem  requires  that  we 
solve  the  known  boundary  problem,  the  numerical  techniques  for  known  problems 
find  use  in  the  unknown  domain.  In  addition  to  using  the  known  techniques,  we  are 
able  to  change  the  unknown  problem  by  transformations  before  their  application  if 
we  so  choose. 

If  we  do  not  manipulate  the  unknown  boundary  problem,  then  we  must  explic¬ 
itly  approximate  the  boundary  throughout  our  calculations.  This  approach  is  called 
the  trial-free-boundary  method  for  free  boundary  problems  and  is  called  the  front¬ 
tracking  method  for  moving  boundary  problems.  We  will  discuss  these  methods  later 
in  this  thesis. 

If  we  do  transform  the  problem  before  trying  to  solve  it,  then  we  temporarily 
remove  the  unknown  boundary.  This  introduces  some  complications,  however.  The 
transformed  problem  will  be  nonlinear  since  the  original  unknown  boundary  problem 
was  nonlinear.  And  we  have  to  perform  a.  recovery  transform  to  get  back  the  unknown 
boundary  at  the  end  of  the  calculations.  The  types  of  transformations  we  could  use 
are  called  front-fixing,  analytical,  and  fixed-domain. 

In  the  front-fixing  method,  the  unknown  domain  is  transformed  onto  a  known 
domain  with  a  more  complex  differential  equation  and  boundary  conditions.  An 
example  of  this  method  is  the  isotherm  migration  technique. 

In  the  analytical  method,  we  use  techniques  such  as  conformal  mapping  to 
obtain  a  new  problem  such  as  an  integral  equation.  Helmholtz  used  conformal  map- 
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ping  to  solve  fluid  flow  free  boundary  problems  in  1868.  Ilis  work  was  extended 
by  Kirchhoff  in  1860  (5:614).  In  1986,  Elcrat  and  Trefethen  made  their  techniques 
the  basis  of  an  efficient  and  effective  method  for  solving  integrals  that  are  like  those 
found  in  Schwarz- Ghristoffel  mappings  (6:251-265). 

The  fixed-domain  method  uses  a  weak  or  generalized  solution  defined  on  a 
known  domain  that  implicitly  has  information  about  the  unknown  boundary.  Ex¬ 
amples  of  this  method  are  variational  inequalities  and  the  enthalpy  method  for  Stefan 
problems.  Baiocchi  used  variational  inequalities  to  study  porous  flow  problems  in 
earthen  dams  (5:618).  The  enthalpy  H  denotes  the  total  heat  content  in  a  region. 
The  enthalpy  method  can  be  applied  to  “mushy”  regions  that  are  the  result  of  the 
mixture  of  two  phases  such  as  water  and  ice.  This  method  is  advocated  by  Solomon, 
Alexiades,  and  Wilson  (20:8—12). 

This  thesis  will  solve  a  moving  boundary  problem  in  its  original  form.  The 
trial-free-boundary  method  and  the  front-tracking  method  often  need  to  be  applied 
to  domains  with  curved  boundaries.  The  approaches  which  are  the  most  popular 
use  integral  equations,  finite-differences  with  boundary-fitted  coordinates,  or  finite 
elements.  The  boundary  element  method  we  will  apply  will  generate  a  Volterra 
integral  equation  which  will  allow  us  to  integrate  around  the  boundary  of  the  domain 
instead  of  a  double  integral  over  the  area  of  the  domain  and  will  thus  simplify  the 
calculations. 

Rubinstein  provided  existence  and  uniqueness  proofs  for  the  one-dimensional 
Stefan  problem  with  general  initial  conditions  in  1947  (18:37-54).  Originally,  lie 
transformed  the  domains  for  each  phase  to  the  interval  [0. 1]  by  reducing  the  problem 
to  a  system  of  nonlinear  integral  equations  of  mixed  type  (Fredholm  with  respect 
to  space  and  Volterra  with  respect  to  time).  This  system  could  then  be  solved  by 
Picard's  method  of  successive  approximation.  However,  it  required  the  evaluat  ion  of 
double  integrals  and  was  therefore  not  effective.  Again,  in  1947.  he  gave  a  method  of 
reducing  the  Stefan  problem  to  integral  equations  of  Volterra  type  based  on  direct 
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use  of  the  heat  potential  (19:217-220).  Existence  and  convergence  was  guaranteed 
in  the  small;  i.e.,  in  some  neigborhood  of  t  =  0. 

Kolodner  (12:1-31)  used  integral  equations  to  solve  the  Stefan  problem  of  the 
freezing  of  a  finite  depth  lake  in  the  1950’s.  He  incorporated  simple  Green's  functions 
in  his  solution.  He  noted  that  the  cases  of  physical  interest  generate  Vollerra  integral 
equations  of  the  second  kind.  These  integral  equations  can  be  solved  by  numerical 
methods  even  though  they  possess  difficulties  near  t  —  0. 

A vner  Friedman  carried  out  detailed  research  in  the  problems  of  the  evapo¬ 
ration  of  a.  spherical  drop  (8:19-66)  and  the  dissolution  of  a  gas  bubble  in  a  liquid 
(9:327-345).  He  also  showed  the  application  of  the  maximum  principle  to  the  Stefan 
problem  (10:201-211),  Later,  he  compiled  a  text  on  parabolic  partial  differential 
equations  including  a  separate  chapter  on  the  Stefan  problem  (7).  He  demonstrates 
how  a.  parabolic  differential  equation  can  be  transformed  into  a  Vol terra  integral 
equation  which  can  then  be  solved  by  the  method  successive  approximations.  This 
is  the  approach  that  will  be  used  in  this  thesis.  In  addition,  we  will  utilize  the 
boundary  element  method  in  performing  the  integrations  required  for  the  solution 
of  the  Stefan  problem. 

Stefan  problems  are  essentially  nonlinear  because  of  the  condition  on  the  mov¬ 
ing  boundary,  but  linearity  usually  exists  in  each  of  the  domains  on  both  sides  of  the 
boundary  and  we  can  still  utilize  the  integral  equation  approach.  It  has  the  advan¬ 
tage  that  only  the  values  of  the  unknowns  on  the  boundary  enter  into  the  solution 
when  it  is  coupled  with  the  boundary  element  technique. 
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II.  Mathematical  Development 


2.1  The  Heat  Equation 

The  partial  differential  equation  which  we  are  trying  to  solve  is 

cpji  =  v'  (A’°V“) + Q 


where  c  is  the  specific  heat  of  the  material,  p  is  the  mass  density,  u  is  the  temperature, 
t  is  time,  A'o  is  the  thermal  conductivity,  and  Q  is  the  heat  energy  generated  within 
the  domain  of  interest.  If  there  are  no  sources  or  sinks  within  the  domain  then 
Q  =  0.  Also,  if  we  assume  a  homogeneous  isotropic  material  with  constant  thermal 
properties  then  the  constants  can  be  combined  into  a  single  constant  I\  =  ho/cp. 
We  arrive  at  the  heat  equation 


For  one  dimension  the  heat  equation  becomes 


d u  _  d2u 

Tt~hd7> 


And  for  two  dimensions  we  have 

d u  _  r.d2u 

m=h  W  + 

In  order  to  solve  the  heat  equation  on  a  domain,  we  must  supply  initial  and  boundary 
conditions. 

2.2  Boundary  Conditions 

In  solving  differential  equations  it  is  not  sufficient  to  state  only  the  equation 
itself.  A  particular  equation  could  have  an  infinite  number  of  solutions.  In  order  to 
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select  a  single  unique  solution  out  of  this  infinite  number,  it  is  necessary  to  specify 
initial  and/or  boundary  conditions  for  the  equation.  Attempting  to  fit  a  solution  to 
the  conditions  can  be  as  difficult  as  trying  to  solve  the  differential  equation  in  the 
first  place. 

Initial  conditions  are  those  that  must  1°  atisfied  by  the  solution  throughout 
the  domain  at  the  instant  when  consideration  of  the  system  begins.  A  typical  ini¬ 
tial  condition  wifi  prescribe  both  the  solution  at  the  beginning  time  and  the  time 
derivatives  up  through  order  m  -  1  at  the  beginning  time.  Here  m  is  the  order  of 
the  highest  time  derivative  in  the  differential  equation. 

Boundary  conditions  specify  the  value  of  the  solution  on  the  boundary  of  the 
domain  (Dirichlet  condition),  the  normal  derivative  of  the  solution  on  the  boundary 
of  the  domain  (Neumann  condition),  or  a  combination  of  value  and  normal  derivative 
on  the  boundary  (Robin  or  mixed  condition). 

The  prescribed  initial  and  boundary  conditions,  together  with  the  coefficients 
and  any  inhomogeneous  terms  in  the  partial  differential  equation,  comprise  the 
“data"  in  the  problem  modeled.  The  solution  depends  continuously  on  the  data 
if  small  changes  in  the  data  produce  correspondingly  small  changes  in  the  solution. 

A  problem  that  is  modeled  by  the  partial  differential  equation  is  said  to  be 
“well-posed"  if: 

1.  A  solution  to  the  problem  exists. 

2.  The  solution  is  unique. 

3.  The  solution  depends  continuously  on  the  data. 

If  any  of  these  conditions  is  not-  satisfied,  then  the  problem  is  said  to  be  “ill-posed." 

Because  the  heat  equation  is  evolutionary,  it  is  necessary  to  specify  an  initial 
condition.  However,  it  involves  only  the  first  time  derivative  of  the  temperature. 
Therefore,  it  will  require  specifying  only  the  value  of  the  solution  at  the  initial  lime 


value  and  not  the  values  of  an}’  time  derivatives.  If  we  specify  the  value  of  the 
solution  at  any  time  other  than  the  initial  time,  the  problem  may  not  be  well-posed. 
This  is  related  to  the  fact  that  it  is  difficult  to  solve  the  heat  equation  backwards  in 
time.  At  any  time  t  >  t  initial  ■>  the  solution  to  the  initial  value  problem  at  an  arbitrary 
point  in  the  domain  depends  on  “all"  of  the  initial  data.  This  implies  an  infinite 
speed  of  propagation  of  effects. 

The  specification  of  the  initial  condition  along  with  the  behavior  of  the  solution 
at  infinity  constitutes  a  Cauchy  problem  for  the  heat  equation.  A  solution  to  the 
Cauchy  problem  is  infinitely  differentiable  with  respect  to  position  and  time  for  each 
point  in  the  domain  and  for  all  values  of  time  greater  than  the  initial  instant.  This 
shows  the  smoothing  property  of  the  heat  evolution  operator.  A  sectionally  contin¬ 
uous  initial  stale  can  always  evolve  forward  in  time.  However,  if  it  is  not  infinitely 
differentiable  with  respect  to  both  position  and  time,  then  it  cannot  have  originated 
from  an  earlier  state.  Thus,  the  heat  equation  is  irreversible  in  that  “forward"  time 
is  distinguishable  from  “backward”  time.  This  property  in  the  mathematical  model 
corresponds  with  the  second  law  of  thermodynamics  in  the  physical  world  modeled. 

Thus  we  see  that  we  must  have  the  correct  number  of  initial  and  boundary 
conditions  of  the  correct  type  in  order  for  our  problem  to  be  well-posed  and  to  have 
a  physically  realizable  solution. 

2:3  The  Method  of  Successive  Approximations 

In  order  to  describe  the  method  of  successive  approximations,  we  begin  by 
considering  the  space  C&  of  bounded  continuous  functions.  If  .4  C  R’\  let  V  be  the 
sot  of  all  functions  /  :  .4  — >  Rm.  The  zero  function  is  the  member  of  V  that  maps  each 
;r  €  A  to  the  zero  elementof  R”'.  We  define  addition  on  V  by  (/+</)(.<■)  =  /(•<•)+</(. r) 
for  all  x  €  A  and  for  all  /,</  €  V.  Also,  we  define  multiplication  by  a  scalar  on  V 
by  (A/)(.r)  =  A(/(.r))  for  each  A  €  R  and  /  €  V.  Then.  V  is  a  vector  space. 
We  now  consider  the  subspace  C  of  V  defined  by  C  =  {/  €  V  |  /  is  continuous}. 


The  subspace  C  of  continuous  functions  is  also  a  vector  space.  We  now  choose 
the  subspace  Cb  of  C  to  be  the  elements  of  C  that  are  bounded,  i.e.,  Cj,  =  {/  € 
C  |  ||/(:r)||  <  M  for  all  x  G  A ,  M  constant  G  R}.  If  A  is  compact,  then  Ci}  =  C 
since  continuous  mappings  on  a  compact  set  are  bounded.  The  new  space  Cb  is  also 
a.  vector  space.  We  define  a  norm  on  Cb  by  ||/|[  =  sup{|/(.r)|  |  .v  €  /l)  for  all  f  G  C;,. 
This  norm  will  exist  since  /  is  bounded.  Every  Cauchy  sequence  in  Cb  converges 
to  a.  point  in  Cb,  so  Cb  is  complete.  Thus,  Cb  is  a  complete  normed  space  which  is 
equivalent  to  saying  that  it  is  a  Banach  space. 

A  contraction  mapping  is  any  map  T  :  C&(/1,  Rm)  — »  Cj,(.4,Rm )  such  that  B  a 
constant  A  G  R,  0  <  A  <  1  and  \\T(f)  -  T(</)||  <  A||/  -  </||  for  all  }\<j  G  Cb(A,R'n). 
The  mapping  T  is  continuous  and  there  is  a  tiniciue  point  in  Cb  that  is  mapped  to 
itself,  i.e.,  3  /0  G  C(,(/l,R"1)  such  that  T(fo)  =  fo-  The  point  fo  is  called  a  fixed 
point  of  the  mapping.  We  now  form  the  sequence 

;i'o  =  / 

•n  =  T{f) 

=  T(T(/))  =  TV) 

i3  =  nnnm  =  tv) 


This  sequence  is  a  Cauchy  sequence  and  it  converges  to  the  unique  fixed  point  of  the 
mapping  T.  Thus,  we  see  that  if  we  have  a  system  which  is  described  by  a.  contractive 
mapping  or  transformation  that  we  are  guaranteed  to  converge  to  a  unique  solution 
if  we  take  the  output  from  the  system  and  apply  it  again  to  the  input  in  a  repetitive 
feedback  process.  This  process  is  illustrated  in  Figure  (2.1).  This  is  the  essence 
of  the  method  of  successive  approximations.  We  determine  an  initial  “guess”  for 
the  solution  to  our  contractive  system.  We  then  input  this  initial  solution  into  the 
system  to  obtain  a  new  solution.  This  new  solution  is  then  used  as  input  to  the 
system  and  this  process  is  repeated  until  we  have  achieved  the  required  convergence. 


We  can  illustrate  the  method  by  applying  it  to  the  solution  of  an  integral 
equation.  Consider  the  equation 


/(*)  =  «+  [*  (2.D 

Jo 

This  equation  is  a  transformation  on  the  function  /  since  we  can  write  it  as 

Tif  )(*)=«+  [vk(x,mm 

Jo 

The  transformation  T  will  be  contractive  and  there  will  be  a  unique  solution  on  the 
interval  [0,  r]  if 

sup  [  \h{x.^)\d(  =  A  <  1 
x€(0,)-]  JO 

This  places  a  restriction  on  the  kernel  of  the  integral  equation  transformation  fc(.r,£). 
The  solution  will  be  the  function  /(; v)  that  is  the  fixed  point  of  the  transformation. 
We  solve  equation  (2.1)  with  an  appropriate  kernel  by  making  a  trial  guess  for 
/(£)  and  substituting  it  into  the  right  hand  side  of  the  equation.  We  perform  the 
operations  required  by  the  equation  and  obtain  a  new  function  /(. r)  on  the  left  hand 
side  of  equation  (2.1).  This  new  function  can  then  be  substituted  back  in  the  right 
hand  side  and  the  process  repeated.  As  the  process  is  repeated,  we  compare  the  new 
function  created  with  the  old  function  that  was  input.  When  they  are  as  close  as  we 
require,  the  repetitive  cycle  is  stopped  and  we  have  converged  to  our  solution  f(x). 
We  will  show  later  how  to  apply  this  method  to  the  solution  of  the  Stefan  problem 
in  one  and  two  dimensions. 

'2. 4  The  One  Dimensional  Problem 

Our  first  example  is  one  dimensional  diffusion  or  heat  propagation.  The  prob¬ 
lem  considered  in  this  section  is  consistent  with  the  problem  presented  in  (7:210) 
where  existence  and  uniqueness  results  are  presented.  Suppose  F  satisfies 


(Jt  —  Uxx 

u{o,t)  =  m 

u(x,  o)  =  m 


0  <  a:  <  ${t ),  t  >  0 

/(<)>0,  t>0 

<j>(x)  >0,  0  <  .r  <  6, 
<j>(b)  =  0,  b  >  0 


U{s{t),t)  =  0  t  >  0,  s(0)  =  b 

UMt).t)  =  -—«(*)  i>i 


The  boundary  condition  for  the  end  x  =  0  is  non-homogeneous.  The  boundary  con¬ 
dition  on  the  moving  boundary,  x  =  s(t),  is  homogeneous  (Diriehlet).  The  initial 
condition  is  non-homogeneous.  We  know  from  physical  constraints  that  the  temper¬ 
ature  in  the  domain  must  be  bounded.  This  set  of  conditions  along  with  the  partial 
differential  equation  form  a  well-posed  problem.  The  domain  of  interest  for  this 
problem  is  shown  in  Figure  (2.2).  To  solve  this  problem,  we  use  Green’s  Theorem 
as  follows.  We  suppose  that  U  satisfies  the  equation 


Ur  =  C'« 


and,  following  (7:220),  choose 


G(;M;£,t)  =  A'(<M;£,r)  -  A'(-;M,£,r) 


where 

~  (g-Oa‘ 
.  4(i-r). 

so  G{x, /;£,  r)  satisfies  the  adjoint  equation 


GT  =  —6'^ 


A'(;M;£,r)  = 


1 


2tt1/2(<  _  r)1/2 


exp 
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Figure  2.J.  Feedback  iu  a  Contractive  System 


Figure  2.2.  Domain  for  the  One  Dimensional  Problem 
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with  the  boundary  condition  at  £  =  0 


6'(s,t;0,r)  =  0 

Then,  we  introduce  Green's  identity 

0  =  {GVt  -  UG&  -  (GU)r  (2.2) 

and  we  integrate  equation  (2.2)  with  respect  to  £  from  0  to  s(r)  and  with  respect  to 
r  from  0  to  t  —  e  to  obtain 

0  =  Jl"  J>«  ~  UG()(  ~  lGU'i'W,h 

rt-t  f»(r ) 

=  l  l  m-uot)*** 

rt-t  rb 

-  /  ( GU)Td(dT 

Jo  Jo 

rt-t  r*(r) 

-/„  l  {Gl'^dr 

=  J^(GU(-VG()\(..(,)<lr 

-  r\GU(-UG()\(^,<lr  (2.3) 

JO 

rb  rt-t 

-  (GU)TdTdS 

Jo  Jo 

r»(t-t)  rt-t 

-  /  (GU)T  drd( 

Jb  J  *->(*) 

In  equation  (2.3)  we  have  used  the  fact  that  for  any  integrable  function  F(£,t). 

rt-t  rs(r)  rt-t  rb 

/  /  F{^T)d£dT  =  /  f  F(Z.T)d£dT 

Jo  Jo  Jo  Jo 

rt-t  f$(r) 

+  Jo  Jh  F(£,T)d£dr 
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and  have  interchanged  the  order  of  integration  in  the  last  two  integrals  in  (2.3).  We 
are  guaranteed  that  we  can  change  the  order  of  integration  by  Fubini's  Theorem 
which  states: 

/{/ F(t>T)<%}dT  =  j  [j  <lZ 

This  relation  holds  for  Lebesgue  integrals  whenever 

/ / 1 F(t>T)\d{dT  <  co 

Fubini's  Theorem  applies  to  Riemann  integrals  when  they  exist.  If  we  take  into 
account  that  U{s(r),  r)  =  0  and  that  G(x,t;  0,  r)  =  0,  equation  (2.3)  becomes 

0  =  /  (?(*,<;  s(T),r)r^(<(r),r)rfr 

Jo 

+  Jo  G((ar,t;  0,r)/(r)dr 

-  t  G{x,t;t,t-€)V&t-€)<lt 

Jo 

+  [b  G(x^S,omyO)d(  (2.4) 

Jo 

-  fj"'  G(x, t-,(,t-e)U((, (-.)</{ 

Jo 

Combining  the  third  and  fifth  integrals,  taking  into  account  that  U(£,  s-1(0)l£>(,  = 
0,  and  taking  the  limit  as  e  — ►  0  we  obtain 

C(-M)  =  jf  G'(;r,t;s(r),T)t^(s(r),r)r/r 

+  /  G'{(.r,t.0.r)/(r)(/r  (2.5) 

jo 

+  I’ G(jrMMT(WW 

Jo 
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If  we  take  —  of  both  sides  of  equation  (2.5)  and  take  the  limit  as  x  — *  -s(t)  and  let 

Vil/ 


U*.()U  = 


we  obtain 


X=j(f) 


V’0M)U,(<)  =  ^V{v,t)\x=l{t)  +  Jo  Gx{${t),t;s{T),T)V{${T),T)(lT 

+  Jo  G fxds(t),  t]  0,  r )/(r )dr  (2.6) 

+  [bGM*)M-omw 

Jo 

(2.7) 


For  details  concerning  the  limit  we  have  taken  see  Appendix  B.  We  now  integrate 
the  boundary  condition 

l',M <M)  =  -jy(<)  <>° 


from  r  =  0  to  r  =  /  to  obtain 


s {t)  =  b  -  k  /  V"(s(t),  r)</r 
Jo 


(2.8) 


where  we  have  used  the  facts  that  s(0)  =  b,  k  is  a  constant  independent  of  t,  and 


We  use  equations  (2.6)  and  (2.8)  as  the  key  part  of  the  solution  algorithm  for  the 
one  dimensional  Stefan  problem. 
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Algorithm  1: 


1.  Guess  s0{t),  Vo(jc,<)U,(t). 

2.  Substitute  into  the  right  hand  side  of  (2.6)  and  (2.8)  to  obtain  S\ (t),  Vj(.r,  <)|a._,(M 
from  the  left  hand  side  of  (2.6)  and  (2.8). 

3.  Iterate  until  convergence. 

4.  Substitute  sn{t)  and  Vn(t)  into  (2.5)  to  obtain  U(x\t). 

This  algorithm  is  the  embodiment  of  the  method  of  successive  approximations.  In 
step  (1),  we  generate  the  initial  input  for  the  method.  In  step  (2),  we  apply  the 
input  to  the  system  in  order  to  allow  the  contractive  transformations  to  give  us  a 
“better”  estimate  of  the  solution  to  the  system.  We  then  apply  this  better  estimate 
to  the  uystem  input  in  a  feedback  iteration  loop  in  step  (3).  We  repeat  until  the 
output  from  the  system  is  within  the  accuracy  envelope  of  the  input  that  we  have 
determined  we  require.  We  are  guaranteed  that  we  will  obtain  convergence  since 
the  transformations  are  contractive.  We  then  use  the  solutions  we  obtained  for  the 
moving  boundary  position  ${t)  and  the  heat  flux  V(t)  by  the  method  of  successive 
approximations  in  the  equation  for  the  temperature  U{x,t)  in  step  (4). 

2.5  The  Two  Dimensional  Problem 

We  now  consider  the  two  dimensional  problem..  Suppose  that  U  satisfies 


Ut 


V(.v,y,  0) 

dx 

dt 

<ji/ 

dt 


a*(U„  +  U„) 
0 


~hPyU\!J.t) 


(*,y)€ft(0>  0  <t<T 
{x\y)  <E  T(0  =  dU(t) 
(JMf)€fl(0) 

(.r,f/)  €  r(f) 

(x,y)er(t) 


(2.9) 
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We  see  that,  we  have  a  homogeneous  (Dirichlet)  boundary  condition  for  the  temper¬ 
ature  on  the  boundary.  The  initial  condition  is  non-homogeneous.  Again,  we  know 
from  physical  constraints  that  the  temperature  in  the  domain  must  be  bounded. 
We  also  require  that  the  temperature  be  bounded  at  the  origin  x  =  y  =  0.  This 
will  lead  to  our  expressing  the  theoretical  solution  to  the  fixed  boundary  problem  in 
terms  of  a  Bessel  function  of  the  first  kind  of  order  zero.  We  will  give  details  on  this 
solution  in  the  next  chapter.  This  set  of  conditions  along  with  the  partial  differen¬ 
tial  equation  are  a  well-posed  problem.  The  domain  of  interest  for  this  problem  is 
shown  in  Figure  (2.3).  The  velocity  of  the  moving  boundary  is  proportional  to  the 
energy  flow  across  the  boundary  since  the  melting  of  the  ice  is  directly  proportional 
to  the  amount  of  heat  energy  transferred.  Fourier’s  law  states  that  the  heat  flux  is 
proportional  to  the  temperature  gradient: 

$  =  -kVU  (2.10) 

The  minus  sign  is  due  to  the  fact  that  heat  energy  flows  in  the  direction  of  decreasing 
temperature  so  that  the  heat  flow  vector  is  in  the  opposite  direction  to  the  tempera¬ 
ture  gradient.  Therefore,  the  velocity  of  the  moving  boundary  is  proportional  to  the 
temperature  gradient: 


v  =  -kQVU 


—  —  A’o(it4  +  jUy) 


(2.11) 


Since 


^  jlx  jhj 


(2.12) 


we  have: 
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(2.13) 


dx 

dt 

<]y 

dt 


-koUx 
-h  Uy 


To  allow  for  different,  expansion  rates  in  the  f  and  J  directions  we  assume  dif¬ 
ferent  proportionality  constants  h\  and  k-2  and  arrive  at  the  last  two  equations 
in  (2.9).  Thus  these  equations  represent  the  transfer  of  latent  heat  necessary  to 
melt  or  freeze  the  ice.  They  are  known  as  ‘Stefan  conditions’,  These  expressions 
for  dx/dt  and  dy/dt  are  similar  to  expressions  found  in  (16:741-752).  Since  the 
heat  operator  L  —  d/dt  —  «2V2  is  not  self-adjoint,  we  introduce  the  adjoint  operator 
L *  =  —d/dt  —  «2V2  and  obtain  the  adjoint  equation: 


G'r  +  «2(G'^  +  Gm)  —  0 


(2.14) 


If  V  satisfies  (2.9)  and  G  satisfies  (2.14)  then 


0  =  «2[(C(G)(  +  -  (UGfk  -  (UG„), ,]  -  (UG)r 


If  we  integrate  this  equation  with  respect  to  r  from  0  to  t  —  e  and  with  respect,  to  ( 
and  i ]  on  the  region  0(t)  and  let  0(r)  =  f 1(t)  —  0(0),  we  obtain 


0  =  7  {(«2  6']<  +  [a2  UnG}tl}<l(dVdT 

J0  JO(r) 

-  /'"*  [  {(«2  UG&  +  [a2UG,}„  +  (UG)r}dtdv(lT 
Jo  JQ(t) 

=  n  [(i%G^  +  (i2U,,Ghn)d(TdT  (2.15) 

Jo  Jr(T) 

-  /  /  [(tiUG^f^  + (i2UG„)ht)}dcrdT 

Jo  Jl'(T) 

-  n  {UGMZdiidr  -  7  {CGUZdndr 

Jo  JQ(  0)  Jo  JQ{t) 
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In  equation  (2.15)  we  have  used  Green’s  Theorem  §c  A  •  h  ds  =  JRf  V  •  .4  d.rdy 
to  obtain  the  boundary  integrals  on  T(t)  from  the  area  integrals  over  fl(r).  For 
(£,  >/)  €  Cl(t  —  e),  let  f*(£,?/)  be  determined  by  {£,rj)  €  T(t*). 

We  can  now  interchange  the  order  of  integration  in  the  last  two  integrals  in  equa¬ 
tion  (2.15)  and  take  these  two  integrals  to  the  left  hand  side  of  the  equation  to 
obtain 


[  [  \uG)TdTdtdi )  +  I  [  "*  (UG)Tdrd£dn 

JU(Q)Jo 

=  /<_<  /  ((t2G~—  -  a*ir™)  da  dr 

Jo  Jr(r)  \  an  an ) 

In  equation  (2.16),  we  have  used  dU/dn  =  U$n$  +  Unn,r  Thus 

J  {VCUt-,-VG\,*,\dtdv 

+  f.  [(^G')|r=<-«  “  WO)\Tsi*^)]d^<h) 

ft  l  f  rr^\i  i 

=  /  /  a  G—  -  U— )dadr 

Jo  Jr(r)  an  an 


(2.16) 


(2.17) 


In  equation  (2.17),  the  last  term  on  the  left  hand  side  and  the  second  term  on  the 
right  hand  side  vanish  because  U  vanishes  on  T(t).  Thus, 


f  UG\T=t-cd(di;  =  I  G(x, y,  t\ £, y, Q)<f>(^ v)d{dy 
Jn(t-i)  Ju(  o) 

ndU 

a2  G— da  dr  (2.18) 

(r)  an 


We  now  choose  G  to  be 


G'(.r,  f;£,f/,r)  = 


1 


47ra2(t  —  r) 


(2.19) 
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which  is  the  two  dimensional  Green’s  function  for  the  heat  equation  on  an  infinite 
domain.  It  is  the  solution  of  the  heat  equation  with  a  Dirac  delta,  source: 

f)C 

=  «2V2G'  +  6{*  -  my  -  -  r)  (2.20) 

The  Green’s  function  G  is  a  solution  of  the  adjoint  equation  (2.14)  and  if  we  take 
the  limit  as  e  — ►  0  we  obtain  from  equation  (2.18) 


U{x<y,t)  =  f  G'(;r,i/,t;^//,0)^,i/)d^/i/ 

Jo(  0) 


+ 


/  J  a'Gl 

Jo  Jr(r)  on 


du 


dcrdr 


(2.21) 


In  equation  (2.21 )  we  have  included  a  £,  y  subscript  on  the  normal  derivative  D/On  be¬ 
cause  we  are  now  going  to  take  the  partial  derivative  of  both  sides  of  equation  (2.21 ) 

ov 

with  respect  to  the  normal  at  (x,y).  We  take  — - ,  take  the  limit  as  (a* ,y)  -*  !'(/), 

Dnr,y 

DU 

and  define  V(x,y,1)  on  T{t)  by  V{x,y  J)|r(<)  =  -p  then 

On  r{<) 


0 


V(a-,i/,0|r(o  =  J  C  V. 

DG_ 

Gy 


+\v(*,y,t)\rm  +  J/  Jr^/±V<bjT  (2.22) 


Note  that,  for  the  integral  on  T(r),  G  and  V  are  evaluated  at  (a:,  y ,  £(<r,  r),  i/(cr,  r),  r ). 
Thus 


V(*,».<)|no=2  /  ,  TT-GM-Wiv 

''W(O)  C/ll  »•»< 


0(0)  ^nx# 


+2«2  f  j 

Jo  Jl'tr 


OG 

(r)  0nXy 


Vela  dr 


(2.23) 
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By  integrating  the  last  two  equations  from  equations  (2.9)  we  obtain  the  following 
two  equations  which  describe  how  the  boundary  moves  with  time: 


M*)lr(o  =  M0)|r(o)  -  h  [  #,(*, y,  r)\{jCa<y^r(T)dT  (2.24) 

J  0 

=  2/<r(0)|r(o)  -  A-2  /  ty.r,*/,T)|(*^)€r(r)dr  (2.25) 

Here  <r  is  an  index  that  denotes  the  particular  point  on  the  boundary  that  is  of 
interest.  We  use  equation  (2.24)  to  determine  the  x-coordinates  of  the  trajectory  of 
this  boundary  point  designated  by  a  and  we  use  equation  (2.25)  to  determine  the 
y-coordinates. 

In  order  to  solve  for  the  temperature  distribution  on  the  domain  and  the  po¬ 
sition  of  the  moving  boundary,  we  will  employ  the  following  algorithm: 

Algorithm  2: 

1.  Cluess  Vo(;r,|/,f)|r,  a'0(;r,y,0|r,  !to(*»y.<)|r. 

2.  Use  Vo(;P,iM)|r  *n  the  right  hand  side  of  equation  (2.23)  to  obtain  Vj(;r,</,/)|r 
from  the  left  hand  side  of  equation  (2.23).  Use  Vj|p  in  (2.21)  to  calculate 
f(r|r»Uj,|r.  Use  these  values  in  (2.24)  and  (2.25)  to  obtain  aq|r,  i/i|r- 

3.  Iterate  until  convergence. 

4.  Use  converged  values  in  (2.21)  to  obtain  U. 

In  step  (1),  we  generate  the  initial  input  for  the  method.  In  step  (2),  we  apply 
the  input  to  the  system  in  order  to  give  us  a  “better”  estimate  of  the  solution  to 
the  system.  We  then  apply  this  better  estimate  to  the  system  input  in  a  feedback 
iteration  loop  in  step  (3).  We  repeat  until  the  output  from  the  system  is  within 
the  accuracy  envelope  that  we  require.  We  then  use  the  solutions  we  obtained  for 
the  moving  boundary  position  and  the  heat  flux  in  the  equation  for  the  temperature 
U(.v,yJ)  in  step  (1). 
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Figure  2.3.  Domain  for  the  Two  Dimensional  Problem 
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III.  Results 


3.1  Temperature  and  Moving  Boundary  Calculations 

We  choose  the  domain  at  time  t  =  0, 17(0),  to  be  the  open  disk  ;r2  +  y2  <  0.36 
in  the  x-y  plane.  We  also  choose  our  initial  condition  U(x,y,0)  =  <j>(.i\y)  = 
Mcyfi*  +  yl )  where  Jo  is  the  Bessel  function  of  order  zero.  Here  c  is  a  constant  that 
we  need  to  determine.  If  we  choose  k\  —  It?  =  0  in  equation  (2.9)  and  thereby  let. 
the  boundary  be  fixed,  we  will  have  a  well-posed  two  dimensional  initial-boundary 
value  problem  for  the  heat  equation.  This  problem  can  be  solved  by  the  method  of 
separation  of  variables  to  yield  the  solution 

U(x,y,t)  =  exp(-2c2/)Jo(ci/r2  +  y1) 

for  the  cylindrical!)'  symmetric  fixed  domain  (l(t)  =  0(0),  t  >  0.  The  boundary 
condition  requires  the  temperature  to  be  equal  to  zero  on  the  boundary  of  the  domain 
for  all  values  of  time.  Therefore,  the  above  equation  for  U{x,y,t)  is  a.  solution  if  0.6c 
is  the  smallest  zero  of  the  Bessel  function  J0.  The  smallest  zero  of  J0  is  equal  to 
2.4048256  and  we  thus  determine  that  c  is  equal  to  4.0080426.  If  we  substitute 
this  value  for  c  into  our  expression  for  the  temperature,  we  will  obtain  the  initial 
temperature  surface  shown  in  Figure  (3.1).  If  we  now  allow  time  to  progress  from 
t  =  0  to  t  =  0.225,  we  will  obtain  the  series  of  heat  surfaces  shown  in  Figures  (A.l) 
and  (A. 2).  The  exponential  decay  with  respect  to  time  is  readily  apparent  in  the 
figures.  Also,  note  that  the  boundary  does  not  move  as  time  progresses.  We  now 
look  at  the  above  solution  for  the  fixed  boundary  problem  at  the  points  shown  in 
Table  (3.1)  and  also  at  the  origin  (0.0).  We  allow  t  to  take  on  the  values  in  the 
interval  0  <  t  <  0.3.  The  solution  at  the  given  points  is  shown  by  the  solid  curves 
in  Figure  (3.2).  Curve  A  is  the  solution  C(.i\y,t)  as  a  function  of  time  at  the  origin 
(0,0).  Curve  B  is  U{.i\y,t)  for  (a-,t/)  being  any  of  the  eight  points  of  group  B  in 
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Temperature 


Figure  3.1.  Initial  Condition  Temperature  Surface 


Figure  3.2.  Solution  for  a  Fixed  Boundary 
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Point  Nuaber  Angle (radians) 


X-Coordinate 


1 

2 

3 

4 

5 

6 

7 

8 


0. 

0.7853981 

1.5707963 

2.3561945 

3.1415927 

3.9289910 

4.7123890 

5.4977870 


0.1500000 

0.1060660 

0.0000000 

-0.1060660 

-0.1500000 

-0.1060660 

0.0000000 

0.1060660 


1 

2 

3 

4 

5 

6 

7 

8 


0. 

0.7853981 

1.5707983 

2.3561945 

3.1415927 

3.9269910 

4.7123890 

5.4977870 


0.3000000 

0.2121321 

0.0000000 

-0.2121320 

-0.3000000 

-0.2121320 


•I'A'I'X' 


0.2121320 


1 

2 

3 

4 

5 

6 

7 

8 


0. 

0.7853981 

1.5707983 

2.3561945 

3.1415927 

3.9269910 

4.7123890 

5.4977870 


0.4500000 

0.3181981 

0.0000000 

-0.3181981 

-0.4500000 

-0.3181980 

0.0000000 

0.3181980 


Y-Coordinate 

0. 

0.1060660 

0.1500000 

0.1060660 

-0.0000000 

-0.1060860 

-0.1500000 

-0.1060680 

0. 

0.2121320 

0.3000000 

0.2121320 


ffivXO. 


-0.2121321 

-0.3000000 

-0.2121321 

0. 

0.3181981 

0.4500000 

0.3181981 

-0.0000000 

-0.3181981 

-0.4500000 

-0.3181981 


Table  3.1.  Solution  Evaluation  Points 


Figure  3.3.  Domain  Increments  Used  for  Integration 

Table  (3.1).  The  plot  of  the  solution  for  each  point  repeats  on  the  same  curve  due 
to  the  symmetry  of  the  problem.  Curve  C  is  U(x.yJ)  at.  the  eight  points  shown  in 
group  C  in  the  table.  And  finally,  curve  D  is  at  the  eight  points  shown  in 

group  D.  We  now  use  our  algorithm  to  solve  the  fixed  boundary  problem.  We  evaluate 
the  integral  over  0(0)  in  equation  (2.23)  by  dividing  the  x-axis  into  80  increments. 
This  will  result  in  a  division  of  the  y-axis  which  is  allotted  proportionally  as  we 
move  across  the  domain  along  the  x-axis.  We  used  this  proportional  spacing  to  save 
computation  time.  It  is  unnecessary  to  use  a  very  fine  grid  as  we  are  integrating  over 
y  when  the  length  of  the  y-segment  is  short,  as  occurs  towards  the  outer  edges  of  the 
domain  as  viewed  along  the  x-axis.  This  is  completely  due  to  the  specific  geometry 
we  have  chosen  for  our  domain.  If  we  had  chosen  a  square,  then  the  y-axis  would  be 
broken  up  into  the  same  number  of  increments  as  the  x-axis.  We  break  the  boundary 
into  32  boundary  elements,  and  the  time  axis  into  100  increments.  The  spatial 
division  of  the  domain  is  shown  in  Figure  (3.3).  We  now  iterate  12  times  to  reach 


Figure  3.4.  Fixed  Boundary  Time  Evolution 


convergence  in  the  successive  approximation.  The  squares  in  Figure  (3.2)  represent 
the  solution  at  selected  time  values,  again  at  the  above  selected  spatial  point  s,  using 
our  algorithm.  The  squares  fall  on  or  very  close  to  the  curves  which  represent  the 
theoretical  solution.  We  will  analyze  the  errors  later.  If  we  draw  a  three-dimensional 
graph  showing  the  boundary  time  evolution  for  the  fixed  boundary,  we  will  have 
the  result  shown  in  Figure  (3.4).  The  vertical  lines  show  the  time  trajectories  of 
the  boundary  points  used  in  calculating  the  boundary  integral  around  the  fixed 
boundary.  The  circles  represent  the  boundary  at  each  moment  of  time. 

We  next  use  our  algorithm  to  solve  the  moving  boundary  problem  correspond¬ 
ing  to  k i  =  h'2  =  1.0  in  equation  (2.9).  For  this  moving  boundary  problem,  we 
iterated  11  times  to  obtain  the  solution  indicated  m  Figure  (3.5).  The  curves  in 
Figure  (3.5)  are  the  same  curves  that  appear  in  Figure  (3.2)  and  thus  they  represent 
the  theoretical  solution.  They  are  included  for  reference  to  compare  the  solution  of 
the  moving  boundary  problem  with  the  solution  of  the  theoretical  fixed  boundary 
problem.  The  squares  in  Figure  (3.5)  represent  the  solution  obtained  using  our  al¬ 
gorithm  at  the  same  points  as  the  data  points  in  Table  (3.1).  The  squares  which 
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indicate  our  solution  of  the  moving  boundary  problem  show  a.  shift,  as  the  boundary 
moves  which  is  as  we  would  expect.  The  temperature  is  higher  at  a  specific  point  at 
a.  given  time  since  the  boundary  of  the  domain  is  moving  outward  from  the  center 
of  the  domain.  The  boundary  time  evolution  for  the  moving  boundary  is  shown  in 
Figure  (3.6).  Again,  the  vertical  lines  show  the  time  trajectories  of  the  boundary 
points  used  in  calculating  the  the  boundary  integral  around  the  moving  boundary.  If 
we  look  at  the  motion  of  the  boundary  in  two  dimensions,  we  will  obtain  the  graph 
shown  in  Figure  (3.7).  The  decay  of  the  temperature  surface  with  respect,  to  time 
is  shown  in  the  series  of  graphs  in  Figures  (A. 3)  and  (A.*l).  We  can  see  that  the 
temperature  falls  off  at  a  faster  than  linear  rate  and  that  tie  surface  spreads  out  as 
the  boundary  expands.  The  expansion  rate  is  equal  in  the  x-  and  y-direclions. 

We  now  will  allow  different  expansion  rates  in  the  x-  and  y-directions.  We  set 
A'i  =  3.0  and  h2  =  0.5.  We  obtain  convergence  after  12  iterations.  The  expanding 
boundary  will  be  graphed  as  shown  in  Figure  (3.S).  It  is  readily  apparent  that  the 
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Figure  3.6.  Moving  Boundary  Time  Evolution 


Figure  3.7.  Moving  Boundary  for  A-j  =  1.0  and  A--2  = 


Moving  Soondoy  Plot 


Figure  3.8.  Moving  Boundary  for  k i  =  3.0  and  k-i  =  0.5 

boundary  is  moving  much  more  in  the  x-direction  than  in  the  y-direction.  This  case 
is  no  longer  cyliudrically  symmetric  and  hence  represents  a  (rue  two-dimensional 
solution  of  the  heat  equation.  The  decay  of  the  temperature  surface  is  shown  in 
Figures  (A. 5)  and  (A.6).  We  can  see  the  greater  broadening  of  the  surface  in  the 
x-direction  than  in  the  y-direction. 

For  added  emphasis  of  the  different  expansion  rates,  we  set  k\  =  5.0  and 
k-i  =  0.5.  The  two-dimensional  moving  boundary  is  shown  in  Figure  (3.9).  The 
expansion  in  the  x-direction  is  very  pronounced.  The  three-dimensional  boundary 
time  evolution  for  the  moving  boundary  is  shown  in  Figure  (3.10).  The  temperature 
surfaces  are  shown  in  Figures  (A. 7)  and  (A. 8).  Here  the  surface  becomes  much 
t  hinner  in  the  y-direction  as  time  progresses. 

■1.2  lu-eculioit  Times 

The  execution  times  were  determined  as  a  function  of  both  the  division  of  the 
time  axis  and  the  division  of  the  space  axes.  The  division  of  the  boundary  was  allowed 


Moving  Boundary  Plot 


Figure  3.9.  Moving  Boundary  for  Ay  =  5.0  and  Ay  =  0.5 


Figure  3.10.  Moving  Boundary  Time  Evolution  for  Ay  =  5.0  and  Ay  =  0.5 
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to  vary  from  vps  =  22  to  ups  =  32  points  in  increments  of  two  elements.  The  three- 
dimensional  plots  of  the  execution  times  are  shown  in  Figures  (A. 9)  and  (A.iO).The 
execution  times  as  a  function  of  the  division  of  the  time  axis,  nt.  are  shown  in 
Figures  (A.  11)  and  (A. 12).  Here  we  notice  that  the  execution  time  increases  almost 
linearly  with  respect  to  nt.  We  expect  this  since  the  time  division  will  determine 
the  amount  of  time  spent  in  the  integration  loop  for  the  line  integral  and  this  is  a 
single  loop  as  opposed  to  the  nested  loop  for  the  spatial  integrations.  The  execution 
times  as  a  function  of  the  division  of  the  space  axes  is  shown  in  Figures  (A. 13)  and 
(A. 14).  Here  one  must  remember  that  the  division  of  the  y-axis  is  dependent  upon 
the  division  of  the  x-axis.  Therefore,  the  label  “nx”  actually  represents  the  compound 
e fleet  of  the  division  of  both  space  axes.  This  is  why  the  growth  of  the  execution 
time  curve  is  geometric  with  respect  to  nx.  The  spatial  integrations  are  performed 
by  a  nested  loop  which  effectively  represents  the  product  of  the  spatial  dimensions, 
thereby  reflecting  a.  quadratic  growth  rate.  Thus  we  see  that  the  execution  times 
are  almost  independent  of  the  number  of  boundary  elements  and  that  they  depend 
most  strongly  on  the  spatial  division  of  the  domain.  Therefore,  in  order  to  keep 
computational  time  to  a  minimum,  it  is  advisable  to  keep  nx  as  small  as  possible 
and  still  retain  the  accuracy  that  we  require.  We  can  allow  np$  to  be  larger  since  it 
has  a  smaller  effect  on  the  execution  time.  The  value  of  the  time  division,  nt,  will 
determine  the  execution  time  at  an  intermediate  rate.  We  can  consult  the  graphs  and 
pick  values  that  give  reasonable  execution  times  for  the  accuracy  that  we  require. 
The  execution  times  given  are  all  for  the  ELXSI  computer  of  the  Air  Force  Institute 
of  Technology. 
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IV.  Conclusions  and  Recommendations 


4.1  Satisfaction  of  Goals 

We  have  presented  the  mathematical  derivation  of  models  for  phase  change 
problems  in  one  and  two  dimensions  using  the  boundary  element  technique.  We 
then  applied  this  boundary  element  technique  to  a  specific  series  of  two  dimensional 
melting  ice  problems.  We  note  that  our  models  are  also  applicable  to  the  solution 
of  the  diffusion  equation.  This  thesis  has  amply  demonstrated  the  facility  of  using 
computer  graphics  as  a  tool  to  display  numerical  results.  The  sheer  volume  of  (  lie 
results  from  the  calculations  would  overwhelm  our  ability  to  see  the  patterns  of 
meaning  without  the  use  of  computer  graphics. 

4.2  Future  Work- 

This  research  is  ripe  for  further  development.  Viable  areas  for  further  work 
include: 

•  Consider  other  two  dimensional  domains  on  which  to  solve  the  heat  equation. 
These  could  include: 

1.  Convex  and  non-con  vex  shapes. 

2.  Symmetric  and  non-symmetric  shapes. 

•  Extend  the  problem  to  three  dimensions.  For  a  development  of  the  three  dimen¬ 
sional  equations  that  is  similar  to  our  development  in  this  thesis,  see  ( 1 J :  105— 
110). 

•  Allow  the  medium  to  become  non-homogeneous  or  anisotropic.  This  implies 
allowing  A'0,  0.  and  p  to  be  spatially-dependent. 

•  Change  the  boundary  conditions.  Possible  changes  include: 
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1.  Allow  the  temperature  to  be  non-zero  along  parts  of  the  boundary  (non- 
homogeneous  problem). 

2.  Insulate  parts  of  the  boundary  (Neumann  condition). 

3.  Allow  mixed  or  Robin  boundary  conditions  along  parts  of  the  boundary. 

•  Make  the  problem  non-liomogeneous  by  allowing  sources  and/or  sinks  of  heat 
energy  within  the  domain. 

•  Allow  the  medium  properties  Ao,  c,  and  p  to  be  time-dependent. 

•  Utilize  alternate  methods,  such  as  Gaussian  quadrature,  for  performing  the 
numerical  integrations. 

•  Use  other  numerical  methods,  such  as  finite  differences  or  finite  elements,  to 
solve  the  Stefan  problem. 

•  Perform  an  execution  time  comparison  and  analysis  between  the  solutions  using 
alternate  numerical  methods. 

•  Perform  an  error  comparison  and  analysis  between  the  solutions  using  alternate 
numerical  methods. 

•  Establish  the  contractive  nature  of  the  mapping  for  the  two  and  three  dimen¬ 
sional  problems. 

4.3  Remarks 

We  would  recommend  that  effort  be  given  to  finding  a  numerically  efficient 
method  for  performing  the  integrations  used  in  solving  the  moving  boundary  prob¬ 
lems.  This  thesis  required  an  enormous  amount  of  computer  time.  Further  thought 
must  be  given  to  improving  the  execution  speed  and  the  memory  requirements  of 
the  code,  especially  if  more  intricate  geometries  are  to  be  considered.  If  fast  comput¬ 
ing  resources  are  not  available,  the  moving  boundary  problem  becomes  a  formidable 
computational  task.  Also,  the  data  files  generated  are  rather  large  and  it  is  not 


possible  to  keep  many  of  them  oil  disk  storage.  Tape  storage  is  one  possibility.  Also, 
further  thought  could  be  given  to  reducing  the  amount  of  required  information  to 
a.  minimum.  We  kept  only  recent  results,  but,  found  later  that  it  would  be  advan¬ 
tageous  to  be  able  to  refer  back  to  earlier  data.  The  computer  graphics  files  also 
require  a  large  amount  of  storage  space. 

Overall,  this  problem  is  quite  intriguing  and  also  satisfying.  It  is  possible  to 
obtain  good  results  that  fit  well  with  one’s  intuitive  expectations.  Our  appetite  has 
certainly  been  whetted,  and  we  would  like  to  consider  future  work  in  this  interesting 
area. 
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Appendix  A.  Figures 
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Figure  A.l.  Temperature  Surfaces  for  Exact  Solution  (1  of  2) 
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Figure  A.2.  Temperature  Surfaces  for  Exact  Solution  (2  of  2) 


Figure  A. 3.  Temperature  Surfaces  for  h'i  =  1.0  and  k-i  —  1.0  (1  ol  2) 


A-3 


HEAT  SURFACE  AT  TIME  «  0.0375 


HEAT  SURFACE  AT  TIME  =  0.0750 


Figure  A.4.  Temperature  Surfaces  for  k\  =  1.0  and  k?  =  1.0  (2  of  2) 
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HEAT  SURFACE  AT  TIME  =  0.2250 


Figure  A.6.  Temperature  Surfaces  for  ki  =  3.0  and  k-2  =  0.5  (2  oi  2) 


Figure  A.7.  Temperature  Surfaces  for  k\  =  5.0  and  k?  -  0.5  (1  ot  2) 
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HEAT  SURFACE  AT  TIME  -  0.0375 


HEAT  SURFACE  AT  TIME  =  0.0750 


Figure  A. 8.  Temperature  Surfaces  for  h‘i  =  i3.0  and  h>  =  0.5  (2  of  2) 
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Figure  A. 12.  Execution  times  as  a  function  of  nl  (2  ol  2) 


A-10 


Figure  A. 14.  Execution  times  as  a  function  of  nx  (2  of  2) 
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Appendix  B.  Derivation  of  the  Limit 
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We  want  to  show,  for  v(t)  (0  <  i  <  a)  a-  continuous  function,  and  the  moving 
boundary  s(t)  (0  <  l  <  cr)  satisfying  a  Lipschitz  condition,  that  for  every  0  <  t  <  a 


Jis)-.  it  Si  mu, u* = i,(r) 


dr 


Here  we  have  used  the  fundamental  solution  to  the  heat  equation  defined  by: 


2xl/2(t  —  r)1/2 


4(t-r) 


In  the  following  C';,  where  i  is  an  integer,  denotes  a  constant.  We  start  by  defining 
the  integral  /: 

'  ■  L  fiir?)  - K(x' i'  a{T)'T),,T  -  L 


We  can  write  I  as  the  sum  of  two  other  integrals,  l\  and  I 2,  where 

f*  ‘V  ~  *(f)  T.' 


/l  =  L|r^A'(a:’f;5(r)’T)f/r 


and 


h  =  I  -11:  ~  f-  [A'(.r, t\ s(t ), r )  -  /v(j(<),<;  j(r),r)](/r 


f'  -g(0  -  3(t) 

f-5  2(t  -  t) 

The  Lipschitz  condition  on  the  boundary  s(t)  gives  us 


|s(*)~s(r)|  <  C\\t  —  T 


This  implies  that 


V^c'Lw^  =  2C'^ 
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so  \I-i\  will  go  to  zero  as  8  goes  to  zero.  We  now  want  to  look  at  |/i|  and  in  order  to 
do  so  we  introduce  the  integral  Jy  defined  by 


■h  =  J  iff-  \KMs[tMdT 

Jt-S  2(t  —  T) 


If  we  make  the  substitution  z  =  (t  —  t)/(x  —  s(t))2  in  Jy  and  note  that  x  —  s(t)  <  0, 
we  will  have 


J,  =  - 


1 


47T1/2 


1  Jo 


Here  (  =  8/(x  —  s(<))2.  When  x  — >  s(f)  then  (  — >  oo.  Evaluating  the  integral  for  Jy 
we  see  that  Jy  — +  We  look  at  J\  —  I\ 


J'-h~Lk 


1  -e 


(;1‘  —  s(r))2  —  (;r  —  s(/))2 

4«  -  r) 


dr 


If  we  expand  the  numerator  of  the  exponent  and  combine  terms  we  see  that  the 
exponent  is  bounded  by 


4(/~r)  M*)  “  *(r)l(l*  ”  a(0l  +  I'2’  ~  s{t)\)  <  C2(|;c  -  s(t) |  -f  | $(t)  -  s(r)|) 

We  have  used  the  Lipschitz  condition  on  the  boundary  to  obtain  the  last  inequality. 
Since  we  are  going  to  let  8  become  close  to  zero  and  hence  have  r  become  close  to  t. 
and  we  want  x  to  approach  s(t),  we  ^an  assume  that  the  last  inequality  is  less  than 
unify.  This  implies  that  the  expression  in  the  brackets  in  J\  — Iy  will  be  bounded  by 


C'3(|.t  —  s(t)\  +  \s(i)  —  ■s(r)|) 


If  we  use  this  fact  in  Jy  —  Iy  and  use  the  known  inequality  ye  y  <  C.y  for  y  >  0.  we 
will  have 
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Since  we  know  that  J\  — ►  —  -  when  x  — >  s(i),  that  |/^|  <  C$81/2,  and  that  I  =  I1  +  I2, 
we  see  that 

limsup/  +  -  <C$81/2  (B.l) 

j_i(o-o  2 

We  can  also  see  that  |/i|  <  C7.  We  now  define  the  integral  I\0 

We  see  that 

A'o  <  \h\  +  [_s  A'(.t,  *(r ),  r )dr 

The  last  integral  is  bounded  because  of  the  Lipschitz  condition  so  we  have 

A'o  <  C8 


We  now  introduce  the  integral  A'j 


A'i  =  /  ,  ^  ^^A(s(0, <; s(t), r)dr 

J  t  — 


2«-r) 

This  integral  is  bounded  by  a  constant  because  of  the  Lipschitz  condition  on  the 
boundary  $(t).  We  now  look  at  the  equation 


L  =  Jo  V^W^t lK(x’t’s(T)'T)dT  ~  Jq  v(T)-~l I  17^A(s(*),  t] s{t), t)(It 


We  break  L  into  the  sum  of  two  integrals  L  =  L\+  L?  where 


Jt-b  2(i  —  t  )  2(/  —  r) 


L-2  =  jf  t>(r)l .  .,.£^1  Jv (,r, s(r), r)(/r  -  jf  {,(7)^!— ^A'(s(/). /;  .s(r).  r)(/r 
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We  see  immediately  that 


lini  L2  =  0 
*-«(«) 

If  we  write  v{t)  =  v(t)  +  (v{t)  —  v(t))  in  our  expression  for  L\  and  use  (B.l)  and  the 
fact  that  A'0  and  K\  are  bounded  by  constants,  we  see  that 

lira  sup  |Ia  +  \v{t)\  <  C9(61/2  +  sup  \v(t)  -  i>(r)| 

jc— #(()-0  —  t— <^<r<t 

We  thus  have 


limsup|(£i  +  Li)  +  <  C'o(61/2  +  sup  |i>(/)  -  v(r)| 

x— «(()-<>  2  t-t<r<l 

Since  the  left  side  is  independent  of  6  and  the  right  side  goes  to  zero  when  6  — >  0, 
we  have 

limsup|(Ai  +  i2)  +  iy(<)  |  =  0 

JC  — ■»(*)— 0  2 

This  is  equivalent  to 


This  completes  our  derivation.  We  have  followed  closely  the  proof  given  in  (7:217-219). 
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